In [1]:
from __future__ import print_function, division
In [2]:
# This changes the current directory to the base saga directory - make sure to run this first!
# This is necessary to be able to import the py files and use the right directories,
# while keeping all the notebooks in their own directory.
import os
import sys
from time import time
if 'saga_base_dir' not in locals():
saga_base_dir = os.path.abspath('..')
if saga_base_dir not in sys.path:
os.chdir(saga_base_dir)
In [3]:
import hosts
import targeting
import numpy as np
from scipy import interpolate
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy import table
from astropy.table import Table
from astropy.io import fits
from astropy.utils.console import ProgressBar
from collections import Counter
In [4]:
%matplotlib inline
from matplotlib import style, pyplot as plt
plt.style.use('seaborn-deep')
plt.rcParams['image.cmap'] = 'viridis'
plt.rcParams['image.origin'] = 'lower'
plt.rcParams['figure.figsize'] = (14, 8)
plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 16
plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 14
In [5]:
from IPython import display
from decals import make_cutout_comparison_table, fluxivar_to_mag_magerr, compute_sb, DECALS_AP_SIZES, band_to_idx, subselect_aperture
Parts adapted from DECALS low-SB_completeness AnaK overlap.ipynb
In [6]:
hsts = hosts.get_saga_hosts_from_google(clientsecretjsonorfn='client_secrets.json', useobservingsummary=False)
anak = [h for h in hsts if h.name=='AnaK']
assert len(anak)==1
anak = anak[0]
In [7]:
bricknames = []
with open('decals_dr3/anakbricks') as f:
for l in f:
l = l.strip()
if l != '':
bricknames.append(l)
print(bricknames)
In [8]:
base_url = 'http://portal.nersc.gov/project/cosmo/data/legacysurvey/dr3/tractor/{first3}/tractor-{brickname}.fits'
for brickname in ProgressBar(bricknames, ipython_widget=True):
url = base_url.format(brickname=brickname, first3=brickname[:3])
target = os.path.join('decals_dr3/catalogs/', url.split('/')[-1])
if not os.path.isfile(target):
!wget $url -O $target
else:
print(target, 'already exists, not downloading')
In [9]:
bricks = Table.read('decals_dr3/survey-bricks.fits.gz')
bricksdr3 = Table.read('decals_dr3/survey-bricks-dr3.fits.gz')
In [10]:
catalog_fns = ['decals_dr3/catalogs/tractor-{}.fits'.format(bnm) for bnm in bricknames]
decals_catalogs = [Table.read(fn) for fn in catalog_fns]
dcatall = table.vstack(decals_catalogs, metadata_conflicts='silent')
In [11]:
sdss_catalog = Table.read('catalogs/base_sql_nsa{}.fits.gz'.format(anak.nsaid))
In [176]:
for dcat in [dcatall]:
for magnm, idx in zip('grz', [1, 2, 4]):
mag, mag_err = fluxivar_to_mag_magerr(dcat['decam_flux'][:, idx], dcat['decam_flux_ivar'][:, idx])
dcat[magnm] = mag
dcat[magnm + '_err'] = mag_err
dcat['sb_r_0.5'] = compute_sb(0.5*u.arcsec, dcat['decam_apflux'][:, 2, :])
dcat['sb_r_0.75'] = compute_sb(0.75*u.arcsec, dcat['decam_apflux'][:, 2, :])
dcat['sb_r_1'] = compute_sb(1.0*u.arcsec, dcat['decam_apflux'][:, 2, :])
dcat['sb_r_2'] = compute_sb(2.0*u.arcsec, dcat['decam_apflux'][:, 2, :])
In [12]:
DECALS_AP_SIZES
Out[12]:
In [287]:
apmag, apmagerr = fluxivar_to_mag_magerr(dcatall['decam_apflux'], dcatall['decam_apflux_ivar'])
apmagres, _ = fluxivar_to_mag_magerr(dcatall['decam_apflux_resid'], dcatall['decam_apflux_ivar'])
In [289]:
apdiff = subselect_aperture(apmagres - apmag, 'r')
apcolor = subselect_aperture(apmag, 'g') - subselect(apmag, 'r')
apmagx = subselect_aperture(apmag, 'r')
good = ~np.isnan(apdiff)&~np.isnan(apcolor)&~np.isnan(apmagx)
plt.scatter(apcolor[good], apmagx[good], c=apdiff[good], cmap='viridis', alpha=.1, lw=0, s=1)
plt.colorbar()
plt.xlim(-2,5)
plt.ylim(26,15)
Out[289]:
In [290]:
for band in 'ugrizy':
reses = subselect_aperture(apmagres, band, None)
print('Band', band)
for ap, res in zip(DECALS_AP_SIZES, reses.T):
print('Aperture', ap,'has', 100*np.sum(np.isfinite(res))/len(res),'% good')
???? Why are so many of the residuals NaN/infs ???
In [291]:
rs = subselect_aperture(apmagres, 'r')
catnotfin = dcatall[~np.isfinite(rs)]
catnotfin['apmagres_allaps'] = subselect_aperture(apmagres, 'r', None)[~np.isfinite(rs)]
make_cutout_comparison_table(catnotfin[np.random.permutation(len(catnotfin))[:10]],
inclres=True, inclmod=True, inclsdss=False, doprint=False,
add_annotation=['apmagres_allaps'])
Out[291]:
Note that this includes only those with r<22 to ensure there's not a flux effect
In [292]:
dmag_of_ap_distr = {}
for ap in DECALS_AP_SIZES:
rs = subselect_aperture(apmag, 'r', ap)
rres = subselect_aperture(apmagres, 'r', ap)
dmag_of_ap_distr[ap] = dmag = rs - rres
plt.hist(dmag[np.isfinite(dmag)&(rs<22*u.mag)], bins=100, histtype='step', label=str(ap), normed=True)
plt.legend(loc=0)
plt.xlabel('r_flux - r_resid')
Out[292]:
In [267]:
ap = 1.0*u.arcsec
rs = subselect_aperture(apmag, 'r', ap)
rres = subselect_aperture(apmagres, 'r', ap)
dmag = rs - rres
perc = 95
p = np.percentile(dmag[np.isfinite(dmag)&(rs<22*u.mag)], perc)
print('nobjs in', perc,'percentile:', np.sum(dmag.value>p), 'cutoff is', p)
msk = np.isfinite(dmag)&(dmag.value>p)&(rs<22*u.mag)
dcatbadres = dcatall[msk]
dcatbadres['dmag'] = dmag[msk]
dcatbadres['r'] = rs[msk]
make_cutout_comparison_table(dcatbadres[:10],
inclres=True, inclmod=True, inclsdss=False, doprint=False,
add_annotation=['dmag', 'r'])
Out[267]:
In [268]:
#cut out the non-overlap region
dsc = SkyCoord(dcatall['ra'], dcatall['dec'], unit=u.deg)
dcutall = dcatall[dsc.separation(anak.coords) < 1*u.deg]
In [269]:
dsc = SkyCoord(dcutall['ra'], dcutall['dec'], unit=u.deg)
ssc = SkyCoord(sdss_catalog['ra'], sdss_catalog['dec'], unit=u.deg)
threshold = 1*u.arcsec
In [270]:
idx, d2d, _ = ssc.match_to_catalog_sky(dsc)
plt.hist(d2d.arcsec, bins=100, range=(0, 3),histtype='step', log=True)
plt.axvline(threshold.to(u.arcsec).value, c='k')
None
In [271]:
dmatchmsk = idx[d2d<threshold]
dmatch = dcutall[dmatchmsk]
smatch = sdss_catalog[d2d<threshold]
In [272]:
idx, d2d, _ = dsc.match_to_catalog_sky(ssc)
dnomatchmsk = d2d>threshold
dnomatch = dcutall[dnomatchmsk]
In [273]:
plt.figure(figsize=(12, 10))
xnm = 'r'
ynm = 'sb_r_0.5'
ap = 1*u.arcsec
apmag, apmagerr = fluxivar_to_mag_magerr(dnomatch['decam_apflux'], dnomatch['decam_apflux_ivar'])
apmagres, _ = fluxivar_to_mag_magerr(dnomatch['decam_apflux_resid'], dnomatch['decam_apflux_ivar'])
rs = subselect_aperture(apmag, xnm, ap)
rres = subselect_aperture(apmagres, xnm, ap)
dmag = rs - rres
dnstar = dnomatch['type']=='PSF '
dnoext = -2.5*np.log10(dnomatch['decam_mw_transmission'][:, 2])
r0 = (dnomatch[xnm] - dnoext)
sb = dnomatch[ynm] - dnoext
plt.scatter(r0[~dnstar], sb[~dnstar],
c=dmag[~dnstar], lw=0, alpha=1, s=3, label='Glx in DECALS, not in SDSS', vmax=0, vmin=-5,
cmap='viridis_r')
plt.colorbar().set_label('r_ap - r_res [{}]'.format(ap))
plt.axvline(20.75, color='k', ls=':')
plt.xlim(17, 23)
plt.ylim(18, 28)
plt.xlabel(r'$r_{0, {\rm DECaLS}}$', fontsize=28)
plt.ylabel(r'$SB_{0.5^{\prime \prime}, {\rm DECaLS}}$', fontsize=28)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.legend(loc='lower right', fontsize=20)
Out[273]:
now inspect those that are in the upper-left of that plot
In [274]:
msk = (r0[~dnstar]<20.75)&(sb[~dnstar]>24)
cat = dnomatch[~dnstar][msk]
cat['dmag'] = dmag[~dnstar][msk]
In [275]:
p = np.percentile(cat['dmag'][np.isfinite(cat['dmag'])], 10)
catlower = cat[cat['dmag']<p]
print(len(catlower))
make_cutout_comparison_table(catlower[np.random.permutation(len(catupper))[:10]],
inclres=True, inclmod=True, inclsdss=False, doprint=False,
add_annotation=['dmag', 'r'])
Out[275]:
In [276]:
p = np.percentile(cat['dmag'][np.isfinite(cat['dmag'])], 90)
catupper = cat[cat['dmag']>p]
print(len(catupper))
make_cutout_comparison_table(catupper[np.random.permutation(len(catupper))[:10]],
inclres=True, inclmod=True, inclsdss=False, doprint=False,
add_annotation=['dmag', 'r'])
Out[276]:
In [277]:
# things we want to identify automatically
disky_things_from_marla = """
mag ra dec unk
19.4785 354.157853072 0.102747198705 false
19.3311 354.255069696 0.619242808225 false
19.1127 354.284237813 0.160859730123 false
18.6914 354.069400831 -0.115904590235 false
19.3392 354.415038652 -0.0645766794736 false
19.4624 354.114525096 0.00532483801292 false
19.1087 354.534841322 0.436958919955 false
19.0242 354.447705125 0.266811924681 false
19.0354 354.136706143 0.691858529943 false
19.0534 354.53888635 0.236989192453 false
19.2452 354.568916976 0.837946572935 false
19.3481 353.473924073 -0.0912764847886 false
19.268 354.011615043 -0.445983276629 false
19.3429 354.408722681 0.904017267882 true
19.1914 354.851596507 0.401264434888 true
19.2529 353.412613626 0.640755638979 false
19.3848 354.878706758 0.364835196656 false
"""
disky_things_from_marla = Table.read([disky_things_from_marla], format='ascii')
In [254]:
sc_marla = SkyCoord(disky_things_from_marla['ra'], disky_things_from_marla['dec'], unit=u.deg)
idx, d2d, _ = sc_marla.match_to_catalog_sky(SkyCoord(dcatall['ra'], dcatall['dec'], unit=u.deg))
np.sum(d2d < 1*u.arcsec)/len(d2d)
Out[254]:
In [314]:
matchcat = dcatall[idx]
apmag, apmagerr = fluxivar_to_mag_magerr(matchcat['decam_apflux'], matchcat['decam_apflux_ivar'])
apmagres, _ = fluxivar_to_mag_magerr(matchcat['decam_apflux_resid'], matchcat['decam_apflux_ivar'])
rs = subselect_aperture(apmag, xnm, None)
rres = subselect_aperture(apmagres, xnm, None)
matchcat['dmag'] = rs - rres
dmagsigs = []
for dmag, dmagdistr in zip(matchcat['dmag'].T, dmag_of_ap_distr.values()):
msk = np.isfinite(dmagdistr)
dmagsigs.append(np.mean(dmagdistr[msk].value)-dmag/np.std(dmagdistr[msk].value))
print(np.mean(dmagdistr[msk].value).shape, dmag.shape, np.std(dmagdistr[msk].value).shape, dmagsigs[-1].shape)
matchcat['dmag_sig'] = np.array(dmagsigs).T
make_cutout_comparison_table(matchcat[np.random.permutation(len(matchcat))[:10]],
inclres=True, inclmod=True, inclsdss=False, doprint=False,
add_annotation=['dmag', 'r', 'dmag_sig'])
Out[314]:
In [ ]: